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Quasi-rigidity means that one builds a theory for assemblies of grains under a slowly changing 
external load by using the deformation of those grains as a small parameter. Is quasi-rigidity a 
complete theory for these granular assemblies? Does it provide unique predictions of the assembly's 
behavior, or must some other process be invoked to decide between several possibilities? We provide 
evidence that quasi-rigidity is a complete theory by showing that two possible sources of indeter- 
minacy do not exist for the case of disk shaped grains. One possible source of indeterminacy arises 
from zero-frequency modes present in the packing. This problem can be solved by considering the 
conditions required to obtain force equilibrium. A second possible source of indeterminacy is the 
necessity to choose the status (sliding or non-sliding) at each contact. We show that only one choice 
is permitted, if contacts slide only when required by Coulomb friction. 

PACS numbers: 45.70.-n, 81.05.Rm, 83.80.Fg 



I. INTRODUCTION 



Historical Overview 



The foundation of many physical theories is the ob- 
servation that a certain physical quantity is "smaU" . In 
practice, this means that the ratio between two different 
quantities with the same units is much less than unity. 
Once a small quantity has been identified, there are two 
ways of proceeding. First of all, that quantity can be 
set to zero, if one wishes to emphasize other aspects of 
the system. This is what is done when presenting the 
harmonic oscillator to students for the first time: one 
usually sets the dissipation to zero, even though its ef- 
fects are quite important. The second possibility is to 
use the small quantity to linearize the equations. This 
is what one does when one suspects that the quantity, 
though small, plays an important role. An example is 
linear stability analysis. 

In the study of granular materials, an obvious choice 
for a small quantity is the distance particles must move 
in order to activate contact forces. This choice is moti- 
vated by the common observation that in a pile of stones 
or marbles, the deformation of the particles due to the 
stresses put on them is not visible to the naked eye. 
Should these deformations be set to zero, or kept as a 
small parameter? 

At the end of the last century, this question was quite 
controversial, as one can see from browsing through a 
conference proceeding from that time 0. The issue at 
hand was the stress distribution at the bottom of a sand 
pile. Several authors j^, JL A 5j proposed theories where 
the grains were assumed to be perfectly rigid. In this 
way, they could circumvent the question of a stress-free 
reference state. However, these theories were criticized 
on many points For example, it was pointed out 
that continuum mechanics could also account for the ob- 



servations .7, §]■ Another objection was that they used 
arbitrary ways to resolve the problem of contact force in- 
determinacy. This problem arises because it is impossible 
to deduce the contact forces in a static granular pack- 
ing from assuming force equilibrium. There is no unique 
solution; instead many solutions are possible j^, lid Il2| . 
This loss of uniqueness occurs because there are more un- 
knowns (contact forces) than equations (vanishing force 
and torque on each particle). 

But the root of all these objections was the realiza- 
tion that rigid particle theories must be radically differ- 
ent from continuum mechanics. In continuum mechanics, 
the counterpart to the particle deformations is the strain. 
Thus neglecting deformations corresponds to eliminat- 
ing the strain. But strain is a fundamental quantity, 
and eliminating it destroys the entire structure of contin- 
uum mechanics. Continuum mechanics can describe the 
macroscopic behavior of granular materials, and it would 
be quite strange if the best microscopic or grain-level the- 
ory were incompatible with it. It was no coincidence that 
most opposition to rigid particle theories came from engi- 
neers, who are more familiar with continuum mechanics 
than most physicists. 

The recently proposed force network ensemble [Tll | can 
be considered as a modern version of the rigid particle 
theories. Instead of adding assumptions to determine the 
forces, all possibilities are considered. The system is thus 
represented by a point in a high-dimensional space 
.12.] . much like in statistical mechanics. One should also 
mention the widely used numerical method of contact 
dynamics |l3l |. which is also based on the assumption of 
rigid particles. 

This paper is concerned with an alternative approach, 
where the particle deformations are a small parameter 
of central importance. We call this approach "quasi- 
rigidity" . This choice is motivated both by the belief that 
such deformations determine the contact forces in phys- 
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ical systems, and by the desire to propose a theory com- 
patible with continuum mechanics, in view of arranging 
a future marriage between microscopic and macroscopic 
theories. 



B. Quasi-rigidity 

Quasi-rigidity was first proposed as the basis for a 
numerical method [H EE IS 113 13 > but has re- 
cently been explored theoretically by a number of authors 
[lilMEllllllililllllilil. Results of this ap- 
proach include a deeper understanding of isostatic pack- 
ings of frictionless particl es ,19, _20,, 22^ , the microscopic 
orig ins of strain 0, |^ l2^7^the stability of packings 
|24l |25| ■ the relation between softening and sliding con- 
tacts |26|| . and jamming 27 1. Work on a corresponding 
numerical method (2lj, |22|, |23 has also continued. 

In quasi- rigid theories, the state of the packing is given 
by the deformation of the particles at each contact. These 
deformations determine the contact forces, which in turn 
govern the motion of the particles. Finally, the parti- 
cle motion gives the change of the deformations. When 
studying the response of a packing to an external load, 
these deformations must be given initial values, analo- 
gous to specifying a reference state in continuum me- 
chanics. 

Do these theories predict a unique evolution of the 
packing? Two possible sources of non-uniqueness have 
been pointed out . First of all, indeterminacy can oc- 
cur if there exists a possible motion that would not mod- 
ify the contact forces. Such motions are called "floppy 
modes" or "mechanisms" . We call this mechanism inde- 
terminacy. Secondly, indeterminacy can arise due to the 
necessity to choose the status (sliding or non-sliding) of 
each contact. This is contact status indeterminacy. 

In this paper, we show that both types of indetermi- 
nacy do not occur. Mechanism instability appears to be 
a problem because force equilibrium is the usual start- 
ing point for quasi-rigid theories. But force equilibrium 
can be considered as a certain limit of Newton's second 
law. When this is done, one can assess the impact of any 
mechanisms that may be present. Contact status inde- 
terminacy can be eliminated by requiring that the con- 
tacts obey Coulomb friction, and letting contacts slide 
only when necessary. This is what is commonly done in 
numerical simulations. A precise definition of 'necessary' 
will be given later in this paper. 

A third possible source of indeterminacy is opening or 
closing contacts. In this situation, one must consider 
transitions between four different statuses(open, non- 
sliding, and two different sliding directions). One would 
like to be sure that such transitions are always uniquely 
determined. Unfortunately, the methods developed in 
this paper are not sufficient to show this, so that this 
possible source of indeterminacy must be investigated in 
the future. 

This paper is organized as follows. Sec. presents 




FIG. 1: Sketch of a biaxial test. An assembly of two- 
dimensional grains is confined by four walls. Forces Fx and 
Fy are applied to the vertical and horizontal walls. 

an overview of the paper, detailing the questions posed 
in the introduction, and sketching the results of the rest 
of the paper. Readers not wishing to savor the details 
may read this section, and then skip directly to Sec. [3 
for a discussion of the results. The stiffness matrix for 
frictional disks is derived in Sec. IIIII and includes a dis- 
cussion of mechanism indeterminacy in Sec. lIIIEl Sec. II VI 
deals with contact status indeterminacy by showing that 
there is a unique way to choose contact statuses. 

II. SYNOPSIS 

A. The Stiffness matrix 

In this paper, we will deal with an assembly of disks, in- 
teracting via Coulomb friction and subjected to a slowly 
changing force. As a concrete example, consider a bi- 
axial box, where a granular sample composed of disks is 
enclosed in a rectangular box of dimensions L^x Ly, with 
forces Fx and Fy exerted on the walls. These forces vary 
slowly with time, and one measures the resulting move- 
ment of the walls. A sketch of the biaxial box is shown 
in Fig.n 

In Sec. IIIII we define the quasi-rigid limit precisely, 
and show that it leads to a piece-wise linear behavior 
of the packing. Thus time can be divided into intervals 
[ti,ti^i] during which the velocities of the particles are 
linearly related to the change in forces: 



where fext represents the external forces {F^ and Fy for 
the biaxial box), v contains the velocities of the particles, 
and k is called the stiffness matrix. 

The motion is only piece-wise linear because the stiff- 
ness matrix k depends on the contact status. Whenever 
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a contact status changes, therefore, k must be modified. 
Therefore, the times {ti} which define the intervals of 
hnearity are the times when one or more contacts change 
status. 



B. Indeterminacy of mechanism 

What happens when the stiffness matrix has a zero 
eigenvalue? In that case, there exists 7^ such that 
kv, = 0. Any multiple of v» can be added to the solu- 
tion of Eq. (Q, and it would still be a solution. Thus it 
would seem that the theory is incapable of determining 
the amplitude of v,. 

But it important to realize that one obtains Eq. 
under the assumption that the external forces change on 
a time scale that is very long compared to the vibrations 
in the granular packing. The appearance of a zero eigen- 
value corresponds to a diverging time scale of vibration, 
and thus the assumptions leading to Eq. (Q) are not met. 
One must use instead Newton's second law, and in this 
case, the amplitude of v, can be determined. If there is 
no interaction between the mechanism and the external 
force (fcxt • V, = 0), then the mechanism is decoupled 
from the other degrees of freedom, and Eq. 0J can still 
be used. 



C. Contact Status Indeterminacy 

Another source of indeterminacy may arise from the 
dependence of k on the contact status. Each contact 
in the packing may be sliding or non-sliding, and each 
choice leads to a different stiffness matrix. If there are 
M contacts, there are 2^ ways to assign contact status, 
and thus 2*^ possible stiffness matrices, and 2*^ different 
solutions to Eq. (^). Which one is correct? As stated in 
the introduction, there is a unique solution if the contacts 
have Coulomb friction, and slide only when necessary. 
Coulomb friction means that the condition 

F = fiF„-\Ft\>0. (2) 

must be obeyed at each contact. Here, F„ and Ft are the 
normal and tangential components of the contact force. 
The constant /i is the Coulomb friction coefRcient. 

When we say that contacts should slide only when nec- 
essary, we mean that they should slide only when they 
would violate Eq. if they did not slide. Since this rule 
places an important role in this paper, we give it a name: 

Principle of minimum sliding: A contact 
slides if, and only if, remaining non-sliding 
would violate Eq. 

To illustrate this principle, let us consider a block 
pushed against a plane with normal force Fn°^*'^ , as shown 
in Fig. [3 A force tangential to the plane Ff'^^^^ is also 



p(ext) 

^(ext) 



FIG. 2: A block pushed against a plane. 

applied. The contact between the block and the plane ex- 

(C) (O 
erts a normal force Fn and a tangential force F^ ' . Let 

these forces be directed opposite to the external ones, so 
-^jiT*^ = Fj^f^ indicates force equilibrium. Let us suppose 
that Fn'^^^^ is fixed, while is slowly increased from 

zero. As long as F^ ''^^^ < fiFn^^\ we have F^'"^ = Fn^^^ 
and Fj'*^'' — F!f^^\ so the block remains in place, and the 
contact is non-sliding. When F^^^*"^ > /i_F'i°'^*' , the block 
must begin to slide, since the contact cannot cancel the 
imposed tangential force without violating Eq. The 
contact is now sliding and F^'^'' = fiFi'^'^*'\ Now sup- 
pose Ft'^^^^ is decreased again, so that i^t^"''*^ < /Lti^i"''*^ 
The block will de-accelerate, and finally stop. When the 
block stops, the principle of minimum sliding says that 
the contact must become non-sliding. Then the block 
will remain in place as long as F^'^^^^ < ^Fn^"^^ ■ Note 
that if we did not apply the principle of minimum slid- 
ing, we would obtain a nonsensical result: the tangential 
contact force would remain constant at ^Fn^'^'^ and start 
to accelerate the block. The frictional forces would thus 
be doing work on the block, which is a violation of the 
second law of thermodynamics. Thus, there is nothing 
artificial about the principle of minimum sliding. It is 
simply makes explicit what is needed to obtain sensible 
results. It also describes what is done in numerical sim- 
ulations. 

Now let us return to the problem of choosing the con- 
tact statuses in a granular packing. It is helpful to con- 
sider the contact forces as points in the {Fn,Ft) plane. 
The set of forces {Fn,Ft) that obey Eq. Q is shown in 
Fig.Oas a shaded region. This set is called the Coulomb 
cone because it forms a cone when plotted this way. 

As the particles move, the contact forces change, and 
thus trace out continuous trajectories in the {Fn,Ft) 
plane. These trajectories must of course remain within 
the Coulomb cone. For contacts in the interior of the 
Coulomb cone, any motion is allowed (for short enough 
times), since there is no danger they will leave the cone. 
Therefore, all such contacts will be non-sliding by the 
principle of minimum sliding. 

Contacts whose forces lie on the boundary of the 
Coulomb cone \Ft\ = ^Fn are called critical contacts, 
and must be handled carefully, since they may leave the 
Coulomb cone. If they are sliding, they will stay on the 
boundary. But as we saw above, we must also allow them 
to leave this surface and enter back into the interior of 
the Coulomb cone. Therefore, each critical contact could 
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FIG. 3: The possible motion of a contact in the Ft) plane. 
The Coulomb cone is shaded. The contact begins at point A 
within the cone. At B, it reaches the surface of the cone, 
becomes sliding, and moves along the cone boundary. At C 
it becomes non-sliding and moves into the cone interior. 



neglecting particle deformations leads to force indeter- 
minacy. We would then have to ask what that physical 
process could be. One possibility is sound waves. As 
we show below, the quasi-static assumption amounts to 
removing 'fast' processes like sound waves. When a con- 
tact changes status, there is probably a "negotiation" 
between the critical contacts, mediated by sound waves, 
that establishes their status. In the quasi-static limit, 
this period of negotiation becomes a single point in time, 
and it is assumed that the principle of minimum sliding 
suffices to determine the new status. Non-uniqueness of 
the choice of contact status means that the details of this 
negotiation must be taken into account. 



D. Uniqueness of the solution 



be sliding or non-sliding. 

It seems that one could determine the status of the 
critical contacts simply by inspecting the particle veloc- 
ities V. These velocities determine the change in con- 
tact forces, and one can easily determine if these changes 
would cause a critical contact to return to the interior 
of the Coulomb cone or not. However, changing a con- 
tact status also changes the matrix k and thus through 
Eq. the velocities v. These new velocities may require 
changing the status of other contacts, provoking another 
re-calculation of v, etc. Thus it seems one must use an 
iterative procedure. One assigns the contact status in a 
certain way, uses Eq. to calculate the particle veloc- 
ities V, and then begins to check if these velocities are 
consistent with the chosen statuses. If an inconsistency 
is found, the status must be changed, and the procedure 
begins again. This procedure must be continued until a 
solution is found. 

In this paper, we are not concerned with this algo- 
rithm, but rather about the uniqueness and existence of 
a solution. Is it always possible to find a solution? Or are 
there many solutions? Note that it is difficult to investi- 
gate these questions numerically. Even though we must 
only deal with the critical contacts, we are often faced 
with situations where all possibilities cannot be investi- 
gated. For example, it is common to have hundreds of 
critical contacts in numerical simulations involving thou- 
sands of particles. This means that the possible ways to 
choose the status cannot even be numbered with 64-bit 
integers. 

Either the non-existence or non-uniqueness would 
bring up hard questions about the quasi-rigid approach. 
If there were sometimes no solutions, the theory could 
not be applied to those situations. On the other hand, if 
the solution were not unique, the stiffness matrix, com- 
bined with the principle of minimum sliding, would not 
be a complete description of the system. Some physical 
process must decide between the different possibilities. 
This unknown process would have been left out of the 
model, leading to indeterminacy in the same way that 



We now sketch the proof that that there is always one, 
and only one choice of contact status that satisfies the 
principle of minimum sliding everywhere in the packing. 
We begin by defining some terms. Let the state of a 
packing be a way of assigning the status to all the crit- 
ical contacts. To each state belongs a corresponding set 
of velocities v, which can be calculated from Eq. A 
state is locally consistent at contact a if the principle of 
minimum sliding is obeyed at that contact, and locally 
inconsistent otherwise. We also refer to the consistency 
of a contact, which means whether or not the principle of 
minimum sliding is obeyed there or not. A state is glob- 
ally consistent if it is locally consistent at all contacts. 
We are thus concerned with the existence and unique- 
ness of the globally consistent state. 

The proof has two premises. First, we assume that all 
possible states lead to a stable packing. The packing is 
stable if 

v^kv > 0, (3) 

for a certain (large) class of relevant vectors v. The sec- 
ond premise is the observation that the left hand side 
of Eq. Q is not modified by the status of the contacts. 
Thus if we consider two different states X and Y for the 
global contact status, one has 

^=k^v^=k%^, (4) 
at 

where k"^ is the stiffness matrix obtained if one chooses 
X , and is obtained by choosing Y. The corresponding 
velocities are v'^ and v^. 

From Eqs. © and Q it is possible to derive a series 
of inequalities, from which one may deduce the following 
theorem: 

Status change theorem: If the status of any 
set of critical contacts changes, the consis- 
tency of at least one of those contacts must 
also change. 
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This statement is sufficient to prove both existence and 
uniqueness. 

To show this, let us cast the theorem into a different 
form. A particular state corresponds to a Mc-bit binary 
number, S £ {0, 1}^^'=, with each bit corresponding to the 
status of a single critical contact, and Mc is the number 
of critical contacts. For concreteness, let us say ^q. = 1 
if contact a is sliding, and = if it is non-sliding. To 
check the consistency of a given state S of contact status, 
we would construct the corresponding stiffness matrix 
Ic^, solve Eq. Q for v^, and then check for consistency 
at each critical contact. The result of this procedure can 
be represented by a second Mc-bit binary number C = 
C{S), where each bit gives the consistency of a critical 
contact, i.e. Ca = 1 if contact a is consistent, and Cq = 
otherwise. Now, the above theorem can be stated as 
follows: Changing any number of bits of S causes at least 
one of the corresponding bits in C to change. This means 
that no two different values of S can lead to the same 
value of C. There are 2*^= possible choices for S, and 2*^"= 
possible values for C, so each possible value of C must 
be associated with a unique value of S. This applies also 
to C = 111 ... 1, corresponding to global consistency. 

This result can be elegantly stated using a more math- 
ematical language. The process of determining the con- 
sistency of a state defines a mapping C of one Mc-hit 
binary number to another: 



C: {0,1} 



s 



C{S). 



(5) 



Now the status change theorem becomes simply: the 
mapping C is bijective. This proves the existence and 
uniqueness of the consistent choice since C is a mapping 
from {0, 1}*^"= onto itself. 



III. THE STIFFNESS MATRIX 

In this section, we present a derivation of Eq. 
where the motion of the particles is related to the change 
in applied force by the stiffness matrix. Several related 
derivations have already been published 0, |2j |2^ . 
The formulation presented here is distinguished from 
these other works in several ways. First, it incorpo- 
rates sliding contacts, which is necessary to consider the 
uniqueness of the globally consistent contact status. The 
second difference is that Eq. is shown to be a cer- 
tain limit of Newton's second law. This is essential to 
resolving the question of mechanism indeterminacy. On 
the other hand, the simplest possible form of the grains 
- disks in two dimensions - is assumed. 

We first describe how interactions between the grains 
are modeled. We then assemble the quantities introduced 
here into vectors and matrices that describe all particles 
in the assembly. We then insert these results into New- 
ton's second law to obtain equations for the motion of the 
grains. Then a limit of these equations is taken, leading 



to Eq. A discussion of mechanism indeterminacy 

completes this section. 



A. Particle interaction model 

We suppose that the grains interact through cohesion- 
less repulsion and Coulomb friction. We adopt the con- 
vention that a positive normal contact force Fn corre- 
sponds to repulsion. Thus the absence of cohesion re- 
quires 

Fn > 0. (6) 
Coulomb friction means that Eq. ^ is 



Recall that 
obeyed. 

When two grains first touch, two springs are created, 
one in the tangential and the other in the normal direc- 
tion. The springs obey Hooke's law so that the normal 
and tangential contact forces Fn, Ft are proportional to 
the spring elongations I?„, Dt- To this restoring force, we 
add a linear damping to model the dissipation of energy: 



Fn = -KnDn " F^K , Ft = -KtDt -TtVt, 



(7) 



where Kn and Kt are the spring constants, r„ and Ft are 
viscous damping coefficients, and Vn and Vt are the the 
normal and tangential relative velocities. Here, £)„ < 
is interpreted as an overlap. 

The springs are stretched by the relative motion of the 
particles, as long this does not violate Eq. ^ or Eq. 10 . 
When the contact is in the interior of the Coulomb cone, 
any motion is possible, so one has 



dPn 

dt 



Vn 



dPt 
dt 



Vt, 



(8) 



where Vn and Vt are just the relative velocities at the 
point of contact: 

Vn = [vi -Vj) ■ nn, 

Vt = {vi - Vj) ■ it + r.uJt + TjUjj, (9) 

where Vi, Ui, and are the velocity, angular velocity 
and radius of particle i, and i and j label the touching 
particles. The vectors fi and t are unit vectors point- 
ing in the normal and tangential directions, respectively. 
Throughout this paper, capital letters indicate quantities 
concerning contacts, and small letters quantities concern- 
ing particles. 

Now let us consider how to handle sliding contacts. It 
is helpful to define 



V ^ fi^Vn + VtSgnDt. 
Note that if the contact is non-sliding, 



dF 



-{Kt + Vt)V. 



(10) 



(11) 
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where F is defined in Eq. (01 • For contacts on the bound- 
ary of the Coulomb cone, we have F = 0. The sign of V 
determines whether such contacts leave or remain within 
the Coulomb cone when made non-sliding. If < 0, the 
contact will move into the interior of the Coulomb cone 
(F > 0). U V > 0, the point would leave the Coulomb 
cone. The principle of minimum sliding can thus be re- 
formulated: 

Principle of minimum sliding: The status 
'sliding' is consistent if, and only if, y > 0, 
where V is defined in Eq. pO(l . 

When a contact slides, Eq. ^ is still valid, but we set 
Ft = and constrain the spring elongations to change so 
that F — 0. This can be accomplished if we use the first 
equation in Eq. JSJl but replace the second with 



dPt 

dt 



Kn 

n sgn Dt 

Kt 



V„. 



(12) 



Once the contact forces are known, the net force / and 
torque r on each particle can be computed: 



a 



(13) 



where the sums are taken over all the contacts that the 
concerned particle makes with its neighbors, and r is its 
radius. 



B. Matrix formulation 

It is useful to consider the proceeding equations in ma- 
trix form. To do so, we must gather the various quanti- 
ties into vectors. To begin with, we can group the force 
and torque exerted on particle i into a vector / . , and the 
contact forces exerted by a contact a into a vector : 




F„ 



Fa,t 



(14) 



It is often convenient to group these vectors together into 
high-dimensional quantities concerning all the particles 
or contacts in the packing: 



fl,x 

n/fi 



fN,x 
fN,y 



(15) 



and 



Fm.ti 
\ FM.t J 



(16) 



Here N is the number of bodies whose motion must be 
considered, and M is the number of contacts between 
these bodies. In these equations, and throughout this 
paper, boldface vectors will denote quantities concern- 
ing all contacts or particles (i.e. vectors in contact or 
particle space), whereas underscores indicate quantities 
associated with a single particle or contact. 
Eq. (|13|l can now be written 



M 



(17) 



a=l 

where c- is a 3 x 2 matrix 



lx»| 



(18) 



This gives the contribution of contact a to the force ex- 
erted on particle i. The symbol Xia is defined as 



Xia 



1 if particle i is first in contact a, 
— 1 if particle i is second in contact a, 
if particle i does not participate in contact a. 

(19) 

If a particle is "first" in contact a, that means that the 
contact exerts a normal force Fn^afia on it. If it is "sec- 
ond" in contact a, a normal force —Fn^ana is exerted on 
it. For each contact between two grains, one element of 
X is 1, and another is —1. x is also called the incidence 
matrix. 

Eq. UTI) holds for each particle {i = 1...N). All of 
these equations can be written compactly using the def- 
initions in Eqs. ((T3|l and lfTC|) : 



f = cF. 



(20) 



The 3A^ X 2M matrix c can be constructed assembling 
an X M array of the c^^ . 

One can consider Eq. as an equation for the un- 
known contact forces F. However, one almost always has 
2M > 37V, meaning that c has at least 2M - 3iV lin- 
early independent null eigenvectors. Therefore, Eq. H2U|) 
does not have a unique solution. This is the force inde- 
terminacy problem discussed elsewhere in the literature 

Iiliailllkj4]. 

We now continue by gathering the other quantities in- 
troduced in Sec. lIII Al into vectors. Eq. {Tj) can be written 



F„ 



-K„D„ 



V or F 



KD - rv, 



(21) 
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Linear Relation 




Equation 


force on particles cx contact forces 


f 


= cF ^ 


contact forces oc spring lengths 


F = 


= -KD iCT 


change in spring length cx relative motion 


D 


= sv iHSJ 


relative motion oc particle motion 


V 





TABLE I: Summary of stiffness matrix derivation as a chain 
of linear relations. The symbol 'oc' is used to mean "is linearly 
related to" . 



where 



K„ = 



Kn 

Kt 



r„ 

Tt 



(22) 



and K and F are 2M x 2M diagonal matrices containing 
the or the on the diagonal. 
Eq. © can be written as 



N 

E 

1=1 



(Ciafvi, or V 



T 



(23) 



where is the transpose of c . Since the dimension of 
V is larger than that of v, Eq. H23() places restrictions on 
V. Not every vector V G K^^^ is allowed, but only those 
vectors in the range of c-^. Physically, this means that not 
every relative motion is possible, but only those that can 
be generated by moving and rotating the particles. The 
dimension of the range of is at most 3A^. There are 
thus 2M — 3N dimensions in M^^^ that are inaccessible. 
These 2M — 3N dimensions are precisely the null space 
of c [13. 

Finally, the relation between D and V in Eqs. ||HJ) and 
(112(1 requires careful treatment, due to the different possi- 
ble contact statuses. Let S be the set of sliding contacts. 
We define a 2 x 2 matrix that depends on the status 
of contact a. li a ^ §, = 1_, if a g S: 



1 



sgn Dt 

Now the relation between v and D can be written: 



dt 



or^ = S(§)V. 



(24) 



(25) 



S(S) is a block diagonal matrix, with the S_i^ on the di- 
agonal. It is a function of §, as indicated. 

Note that sgnDt in Eq. is a constant. In order 
for sgn Dt to change, the contact must cross the Fn axis 
in Fig. 121 This can only happen if the contact passes 
through the interior of the Coulomb cone. In that case, 
the contact would be non-sliding, and Eq. (|24|l would not 
be applied. The exception to this occurs when a contact 
approaches the origin. This brings up the question of 
what happens when a contact opens or closes. We are 
not dealing with that problem in this paper. 

Our derivation of the global stiffness matrix is summa- 
rized in Table. ^ It rests on a chain of linear relations 
that can be established (or assumed) between the various 
quantities. 



C. Equations of motion 

At this point, most derivations of the stiffness matrix 
proceed directly to force equilibrium, and assume that 
the net forces f exerted on each particle are balanced 
by some externally imposed load fext- We take a longer 
route that gives more insight into the situations where 
force equilibrium does not hold. We begin with New- 
ton's second law, which relates the accelerations of the 
particles to the forces exerted on them: 



m— = f 

dt 



(26) 



Here, m is a diagonal matrix containing the masses and 
momenta of inertia of all the grains. We could also write 



m. 



with 



f rrii 
m.j 



(27) 



(28) 



V h/r 



where rrii is the mass of particle i and li is its moment 
of inertia. 

Combining Eq. (jSOJ, (EU and gives 



dv 
' dt 



cKD + 1 



(29) 



This equation can be differentiated once with respect to 
time, and Eqs. H23|) and (I25|l can be used to obtain 

o'^v y J. dv dfext 
m— 7 = cKSc V -KD c Fc— H ; — . 30) 

dt^ dt dt dt ^ ' 

The combination cKSc"^ appears often, so we define k = 
cKSc"^ and write 



m^5- = kv KD c^Fc — 

dt^ dt dt 



dfor. 



dt 



(31) 



This equation gives the full motion, without approxima- 
tion, of the disks. Such an equation is solved numerically 
in the "molecular dynamics" simulation method. On the 
left hand side is the mass times the acceleration (differen- 
tiated by time), and on the right hand side are the forces 
exerted on the particles (also differentiated by time). 



D. Quasi-static balance 

If one makes the quasi-rigid and quasi-static assump- 
tions, then two terms dominate in Eq. 1(311) . The quasi- 
rigid assumption means that the hardness of the particles 
is assumed to be much greater than the confining pres- 
sure, and the quasi-static assumption is that the external 
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forces fext change much more slowly than any timescale 
associated with the contact forces. We describe the con- 
sequences of each of these assumptions below. 

Let us begin with the quasi-rigid assumption. We will 
compare the first two terms on the right hand side of 
Eq. lp?T|l . From Eq. ifTI^ . we see that || will be propor- 
tional to Carrying out this differentiation: 



dn 
'dt 



d 
'dt 



■ tt 



< 



(32) 



Here, i? is a typical particle radius, and V a typical rel- 
ative velocity. Now we can estimate the sizes of the first 
two terms on the right hand side of Eq. I|31|l : 

kv = cKSV - O (KV) , ^KD - O [KV^ 

dt \ R ^ 

where we have used c ^ S ^ 0(1). The quasi- rigid 
assumption implies that the deformations D are much 
smaller than a typical particle radius R, or D/R ^ 1. 
Thus the second term on the right hand side is much 
smaller than the first, and Eq. (|31() becomes 



d^v rp dv 

m — TT = kv c Fc — 

dt^ dt 



dL, 



dt 



(34) 



The term that we have just neglected is called the "geo- 
metric stiffness" |2^. It can again become important in 
some situations. One example will be given in Scc. lIIIEl 
Now let us proceed to the quasi-static assumption. As 
stated above, this means that the external force fext is 
assumed to change much more slowly than any time scale 
associated with the vibrations in the granular assembly. 
To express the separation of time scales, we write 



d 
'dt 



d 

dto 



dti 



(35) 



where the variable to measures a fast time scale, and ti 
a long time scale. In the problem at hand, the fast time 
scale describes vibrations of the packing, and the slow 
time scale is given by the change in the external load 
fext- The presence of e <C 1 before the derivative with 
respect to ti shows that these derivatives are small, as ti 
measures slow changes. 

The assumption of quasi-staticity can be expressed by 
saying that the particle positions x depend on both t^ 
and ti, but the external force f^xt depends only on the 
slow time ti. 



X = xo(io) + xi(ti), f(,xt = foxt(ti). 
After differentiation by time, we have 

V = vq + evi. 
and the 0{1) terms of Eq. H31() are 



m , „ = kvo cFc — — . 

OTq dto 



(36) 
(37) 

(38) 



This equation resembles that of a damped, harmonic os- 
cillator, with three differences. First, it is an equation 
for the velocity, not the position. Second, it is a vec- 
tor equation, not a scalar one. Finally, k and F depend 
on contact status. Nevertheless, it has the same proper- 
ties as a damped, harmonic oscillator. The matrix m is 
positive definite and c-^Fc is non-negative definite, as m 
and F are both diagonal matrices with positive or non- 
negative entries. The stiffness matrix k, however is not 
necessarily positive definite. If there are vectors v 7^ 
such that v-^kv < 0, then k acts like a negative num- 
ber in Eq. H38|l . and Vq grows exponentially on a short 
time scale. This is "motion through an instability" p^ . 
Physically, the external forces push the particles away 
from the positions they must occupy in order to be in 
equilibrium [.27| . In this case, one does not obtain the 
quasi-static balance; rather the packing is unstable and 
is set in motion. 

On the other hand, if v^kv > then k acts like a 
positive number in Eq. (|38|l . and v undergoes damped 
oscillations. In this case vq ^ as ^ oo- One can 
then examine the 0(e) equations assuming vq — 0. (The 
situation when v"^kv = is more complicated, and will 
be discussed in the next section.) Keeping the 0(e) terms 
of Eq. gives 



dir... 



dt 



— kvi. 



(39) 



This is the quasi-static balance, and the same as Eq. 



E. Mechanism indeterminacy 

Now let us consider mechanism indeterminacy. This 
occurs when k has null eigenvalues, i.e., when kv* = 
for v» ^ 0. The amplitude of these motions cannot 
be determined from Eq. H39|l . thus they appear to be 
a source of indeterminacy, just as null eigenvalues of c 
[see Eq. (|20() ] indicate force indeterminacy in rigid parti- 
cle theories. But there is one crucial difference between 
Eq. ||2ni) and Eq. the matrix k in Eq. ^ is a 

square matrix, so if it has null eigenvalue, its range has a 
lower dimension than the left hand side of the equation. 
Thus there are external loads for which Eq. H39|l has no 
solution. On the other hand, the range of c can have the 
same dimension as the right hand side of Eq. (|20l) , even 
when c has many null eigenvalues. 

What happens when Eq. (|39|l has no solution? In that 
case, the quasi-static assumption used to derive Eq. H39|) 
is invalid. This is because a null eigenvalue corresponds 
to a diverging period of vibration in the packing. Thus 
one cannot assume that the force is changing on a time 
scale much greater than the period of vibration. There- 
fore, one must consider Eq. (|34|) without assuming any 
separation of time scales. If one considers only the eigen- 
vector v», Eq. can be integrated once to give 



c?v* 

^~dr 



2^Fcv, 



(40) 
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If the damping can be neglected, then |y» | oc t|foxt|/|m|. 
This is "motion through a mechanism" |24j . The growth 
of the velocity with time is much more gentle than for 
motion through an instability that was discussed in the 
previous section. 

However, Eq. H39|l can have a solution even when k has 
null eigenvalues. In this case, the imposed force does not 
excite the mechanisms. The mechanisms are irrelevant 
to the evolution of the system, and thus the derivation of 
Eq. I|39|) is again valid. Such irrelevant mechanisms are 
quite common. For example, if one constructs the stiff- 
ness matrix of the biaxial box discussed in Sec. Ill Al it is 
convenient to consider the walls as particles. Then, the 
stiffness matrix automatically has three null eigenvalues: 
two associated with the translation of the whole appara- 
tus and one associated with its rotation. As long as the 
forces on opposite walls are equal, these modes will not 
be excited, and they are irrelevant. Another example is 
a "rattler" - a grain which has no contacts. Each one of 
its degrees of freedom yields a zero eigenvalue of k, but 
it can be removed from the system without changing the 
quasi-static behavior. 

Another thing that can happen when a mechanism is 
present is that motion can be stabilized or destabilized 
by the geometric stiffness ^KD in Eq. H31|l . which was 
neglected in the quasi-rigid hmit. When there is a mech- 
anism, the particle displacements are no longer required 
to be small, so this term needs to be considered. An ex- 
ample is when two elliptical particles are pressed together 
[23 ■ If the particles are circular, there is a mechanism: 
the two particles can roll like bearings relative to one an- 
other. If the particles are made elliptical, the neglected 
term ^KD must be considered to predict the behavior. 

Thus, null eigenvalues of k can always be put into one 
of two classes. If Eq. (|39ll has no solution, the null eigen- 
value signals the collapse of the packing, and the quasi- 
static assumption fails, and Eq. 134|l must be used to 
predict the motion of the grains. On the other hand, if 
Eq. (|39|) still has a solution, even though the range of k 
has been reduced in dimension, the null eigenvalue cor- 
responds to a degree of freedom that is irrelevant. Null 
eigenvalues appear to cause indeterminacy only because 
Eq. (|39|) is considered as the most fundamental equation. 
However, Eq. is in fact an approximation to Eq. 



IV. CONTACT STATUS INDETERMINACY 

In this section, we consider contact status indetermi- 
nacy. We will show that only one state does not lead to 
a violation of the principle of minimum sliding at one or 
more contacts. To do so, we must compare the stiffness 
matrices of the different states. We begin by present- 
ing the two hypotheses needed for the proof: first, that 
all possible states are stable, and second that the applied 
load is independent of the state. Then we consider an ex- 
ample where only two contacts are critical, and show that 
the consistent state always exists and is unique. Then 



we consider the general case with an arbitrary number of 
critical contacts. 



A. Conditions needed to siiow uniqueness 

1. The stability condition 

A pac king is stable if the quadratic form Q = v"^kv is 
positive p4. 125ll27l|: 

Q(v,§)=v'^k(§)v>0, (41) 

where § is the set of contacts that are sliding. As we 
saw in Sec. IIII Dl the motion cannot be assumed to be 
quasi-static when Q < 0. 

The quadratic form plays an important role in this 
paper, so we will discuss how it can be calculated. If we 
recall the definition k — cKSc"^, and group factors in a 
suggestive way, we have 

Q(v,§) = [v^c]^KS[c^v]. = V^[KS]V (42) 

The matrix KS is block diagonal, with each block corre- 
sponding to a contact. Thus the Q(v,§) reveals itself to 
be simply a sum over contacts: 

M 

= Ee^'+E^i'^ (43) 

where Qo^^'' is the contribution of contact a if it is non- 
fs) 

sliding, and Qa is its contribution if it is sliding. Using 
Eqs. (|22J) and (|23I), we have 

(v) = K,Vl^ + KtV,% (44) 
Qis)(v) = KnVl„-fiKnVt,c.Vn,c.sgnDt,a.. (45) 

In the following, it is useful to use Eq. (|10|l and replace 
Vn^a with Va- Eq. (|45|l becomes 

Qi'Hv) = Qi'''^(v) - KtVt.,a.Va^sgnDt,o.- (46) 
Now let us define 

F„ =Xtl/t,„sgnA,a, (47) 
so that Eq. H46I) becomes 

Qi'nv) = Or)-F„K.. (48) 
Therefore, the stability condition Eq. 141|) is 

Q(v,§) = O(v,0)-^/'„14 >0. (49) 

Note that Q(v, 0) > 0, because the contribution of each 
contact must be positive. This means the only way to 
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- Status - 




State 


X 


Y § 


all others 


X 


s 


NS S 


NS 


Y 


NS 


S S 


NS 



TABLE 11: The two states X and Y considered in the in- 
dependent load condition. The status of the contacts in the 
sets X, Y, and § are given (S=sliding, NS=non-sliding). All 
contacts not in these sets are non-sliding in both states. 

obtain an unstable packing is for the sliding contacts to 
make large and negative contributions to Q. 

In the following, it will be necessary to compare Q 
for different states. If the sliding contacts present in a 
given state are divided into two disjoint sets Si and §2 
(§1 n§2 ^ 0), then 

Q(v,§iU§2) = 0(v,Si)- ^ F„I4, (50) 

aeS2 

2. The independent load condition 

Let us consider two different states X and Y , each 
with a different set of sliding contacts. Let § be the 
set of contacts sliding in both states, X be set of sliding 
contacts unique to X and Y be those unique to Y (see 
Table nH) . Let v'^ be the velocities in state X and be 
those in Y. Similarly, k-^ = k(SnX) and = k(§n Y). 
If the externally applied force is independent of contact 
status, 

^ = k^v^ = k^v^, (51) 
at 

or both and are the velocities caused by the same 
external forces, but with different stiffness matrices. 
Eq. (|51|l can be rewritten 

k^v^ - k^v^ = cK [S^V^' - S^V^] = 0. (52) 

Now let us multiply this equation from the left by (v^ — 
v^)^: 

[V^ - V^'] ^ K [S^ - V^'] = 0. (53) 
This again is simply a sum over contacts: 

M 

Y.iy^-y-iV [s^^ - s^Yi] = 0. (54) 

a=l 

There will be four types of contributions, corresponding 
to the four columns in Table UTI 

• Contacts which slide in X but not in Y (the set X), 

• Contacts which slide in Y but not in X (the set Y), 

• Contacts which are sliding in both X and Y (the 
set §), 





Contact status 






Consistency 


State 


P 


7 


13 


7 


requirements 


A 


NS 


NS 


1 




/3, 7 not critical 


B 


S 


NS 


1 


1 


7 not critical 


C 


s 


S 




1 


ai = 1 


D 


NS 


S 




0-2 


(Tl = -1, (T2 = 1 



TABLE III: The states considered in Sec. HVBl Contacts (3 
and 7 are critical, so there are four possible states labeled by 

A, B, C, and D. The table gives the status of each critical 
contact in each state (S=sliding, NS=non-sliding). Also given 
are the signs of V. and are known to be positive, be- 
cause it is assumed that /3 becomes critical while the packing 
is in state A, and 7 becomes critical while in state B. The 
other values are deduced from the one contact status change 
theorem given at the end of Sec. 11 V B II o\^ai = ±1 are 
unknown, but used here to show relations between different 
states. 

• Contacts which are non-sliding in both X and Y . 

For the last two classes of contacts, — Sf , so their 
contributions here will be the same as to the quadratic 
form. For contacts a S X the contribution is: 

Kn [Vlo^-Vlof^ sgn Dt^^-rK,V^^^ ) {V^^ - ) . 

(55) 

Defining = — F^ , this quantity can be rewritten 
as 

QrHv^-v^)-Fi^^y/. (56) 

In the same way, the contribution of contacts a € Y is 

Q'r\-''--'') + FrvI. (57) 
Thus Eq. H53|) becomes 

Q(v^ - §) = 5: F^^V^ - Y: FrV^. (58) 

B. Small numbers of sliding contacts 

In preparation for treating the general case, we will 
consider the problem of a packing that may slide at two 
different contacts /? and 7. The four different possible 
states are shown in Table lllll and labeled A, B, C, and D. 
We will use superscripts to indicate quantities belonging 
to each state. For example are the particle velocities 
in state A and S*^ is the status matrix in state C. 

The system starts in state A with no sliding contacts. 
Then contact f3 reaches the boundary of the Coulomb 
cone and becomes sliding, and the packing moves to state 

B. Then contact 7 reaches the boundary, the system 
moves to either state C, where both /? and 7 slide, or to 
state D, where only 7 slides. 

We will consider the questions of existence and unique- 
ness. For example, when /3 becomes sliding, is it guar- 
anteed that B is consistent? If it were inconsistent, then 
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contact (3 should become non-sliding. But if it became 
non-sliding, the system would return to state and (3 
would leave the Coulomb cone. A solution would not 
exist. The question of uniqueness arises when contact 
7 reaches the boundary. If both states C and D were 
consistent, the system could move to either C or D, and 
the solution would not be unique. We will show that the 
consistent state exists and is unique. 

1. One sliding contact 

We first consider the transition from A to B. This 
transition occurs when the contact /3 reaches the bound- 
ary of the Coulomb cone. It starts somewhere within the 
cone, that is with F/3 > [see Eq. (|2Jl]. As the contact 
moves toward the boundary, decreases and then van- 
ishes when P reaches the boundary. Therefore, Eq. I|ll|) 
requires that 

> 0. (59) 

Thus a 1 is given in the top row of the third column of 
Table Iml 

One usually supposes without comment that the state 
B is consistent. But this is not obvious, because all par- 
ticles change their velocities when the state changes. The 
state B will be consistent only if > 0, and no one has 
shown that this must be so. 

To show that B is indeed consistent, let us apply the 
independent load condition to the transition between A 
and B. Setting AT = A, F = S = X = 0, and Y = {/?}, 
Eq. becomes 

Qiv^^v^,9) = -F^^Vf. (60) 

If state A is stable, the quadratic form must be positive, 
leading to 

i^r^/f < 0' (61) 

which gives us some information about the sign of , 
but also unfortunately involves the unknown quantity 
. More information can be obtained by requiring 
state B to be stable: 

Q(v4^v^,{/3})>0, (62) 

and after using Eq. (|50|l 

Q(yA _ 0) „ pAB^yA _ yB-^ > 0^ (g3) 

and finally using Eq. H6U|) : 

F^^Vft < 0, (64) 

Together Eqs. and ((SU show that V^f and have 

the same sign. We already showed that > in 

Eq. ((Snj, thus V/f > as well. Therefore state B is 
compatible, and the solution exists. 



Before proceeding, let us pause to note that the rea- 
soning we have just employed does not depend on state 
A being without sliding contacts. Define S to be the set 
of contacts sliding in state A. If we simply replace the 
empty set in Eq. l(Hn|) and with §, and {/3} with 
{/?} n S in Eq. H62|) . the reasoning remains unchanged. 
Thus we have a general statement: If two states A and 
B differ only in the status of a single contact (3, then 
sgnV^ = sgnVJ^. This state imphes that the consis- 
tency at contact /3 must be different in states A and B. 
If sgn — sgn Vf-f = — 1 (or if /? is not a critical contact) 
then A will be consistent at /3 and B will be inconsistent. 
If sgnV^ = sgnV^ = — 1 then A will be inconsistent 
and B consistent. Thus we have proved a special case of 
the status change theorem: 

One contact status change theorem: If two 
states A and B differ only in the status of 
a single contact /3, then sgnV^ — sgnV^, 
meaning that they differ in the consistency 
at /?. 

Using this theorem, we can now fill in the third and fourth 
columns of Table IIIII States C and D differ only in the 
status of contact /3, so sgnV^ — sgnV^ — a\. States A 

and D differ only in the status of contact 7, so sgn = 
sgnV^^ = (72 • For the same reason, sgnV^^^ = sgnV^, 
and in the next section, we show sgnV^^ — 1. 

2. Two sliding contacts 

Suppose now that the system is in state B, when a 
second contact 7 reaches the boundary of the Coulomb 
condition. Following the same reasoning as at the begin- 
ning of the preceding Sec. lIVBTl we see 

> 0. (65) 

Thus a 1 is given in the second row of the fourth column 
of Table The one contact status change theorem tells 
us we should immediately put a 1 in the third row of the 
same column. 

Now let us consider the uniqueness of the globally con- 
sistent state. When 7 becomes sliding, the system can 
now move to either state C where both /3 and 7 slide, 
or to D, where only 7 slides. We know that the system 
cannot return to A, because > 0, so that state will 
no longer be consistent. 

The one contact status change theorem shows that 
the solution is unique. Only a single contact is differ- 
ent between states C and D, and the theorem states 
that only one of these states can be consistent at /3. If 
sgnV^ = sgnVg^ = (Ji = 1, then contact (3 must slide, 
and C is consistent but not D. On the other hand, if 
(Ti = — 1, then contact (3 must be non-sliding, and D is 
consistent but not C. 

Now let us check the existence of the solution. Suppose 
tJi = 1, meaning the system must move to state C. In 
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order for this state to be compatible, we must also have 
V^' > 0. As shown above, the one contact status change 
theorem guarantees this. 

What happens if cti = —1? In this case, the system 
must move to state D, but is this state compatible? The 
first order state theorem says sgnV^^^ — sgnV^ — (72, 
but this does not help prove compatibility, since a2 is 
unknown. However, we do know that > 0, so let 

us now search for a way to relate to . This can- 
not be done by the one contact status change theorem, 
because B and D differ at two contacts. Therefore, we 
must return to two hypotheses discussed in Sec. IIV Al 
and deduce more information from them. 

If we apply the independent load condition to states B 
and D, we obtain 

g(v^ - v^, 0) = Fj^^V-/ - ^V;^. (66) 

Requiring all four states A, B, C, and D to be stable 
leads to the inequalities 

F^^^f/>Ff^Kf, (67) 

F^''Vp^>F^^^V^^, (68) 

F/^y/ > Ff ^VCf , (69) 

F^^^t^f > Ff ^Kf . (70) 

Now let us suppose that sgn VJ^ sgnV^ (i.e. ai = — 1). 
This means that the left hand side of cither Eq. H67(l or 
(|68|l must be negative. The right hand side of these two 
relations is identical, and must be less than some negative 
number: 

F^^V^" < 0. (71) 

Similar reasoning with Eqs. H69(l and H7()|l leads to 

Ff ^Vf < 0. (72) 

These two inequalities require that and have the 
same sign. Thus we have shown 

sgnf/ ^sgnT//^ ^ sgnVf=sgnVf. (73) 

In the same way, one can begin with the assumption that 
and have opposite signs, and obtain 

sgnVf^sgnl/^ => sgnf/ =sgnV;f . (74) 

Now let us use these results to determine the consis- 
tency of state D. Using the information given in Ta- 
ble iftH Eq. (ESJ becomes 

<Jl^l ^ (72 = 1. (75) 

Recall that fJi = — 1 means that state C is inconsistent, 
and tTi = — 1, (72 = 1 is the condition required for D to 
be consistent. Thus the consistent state always exists. 

But Eqs. H73|l and H74|l also have a much more profound 
implication. When moving from state B to state D, the 



State 


A\A' 


A' 


B\B' 


B' 


§ 


all others 


A 


S 


S 


NS 


NS 


S 


NS 


B 


NS 


NS 


S 


S 


S 


NS 


C 


NS 


S 


NS 


S 


S 


NS 



TABLE IV: States for the proof of the status change theorem. 
Two states A and B are considered. For the proof, a third 
state C is constructed. Contacts in A are sliding in A, but not 
in B, and contacts in B are sliding in B, but not in A. Sets 
A' C A and B' C B are sliding in state C. (Here S=sliding, 
NS=non-sliding). 

sign of V at both /3 and 7 cannot change. Note that this 
reasoning holds even if there are other sliding contacts, 
as long as they remain sliding in all four states A, B, C, 
and D. Thus we have a second special case of the status 
change theorem: 

Two contact status change theorem: If two 
states A and B differ only in the status of 
two contacts f3 and 7, then sgn = sgn Vj^ 

or sgnV^ — sgnV^ , meaning that the con- 
sistency at /3 or 7 (or both) must change. 

C. Many sliding contacts 

In this section, we prove the status change theorem. 
Let us restate it in this way: 

Status change theorem: If two states A and 
B differ only in the status of n contacts 
/?!,.../?„, then for at least one contact f3i, 
we have sgn Vg^ — sgn , implying that the 
consistency at (3i must also change. 

Recall that in Sec. Ill Dl we showed that this statement is 
sufficient to prove that the globally consistent state exists 
and is unique. We now establish this theorem. 

Consider two different states A and B. The contacts 
in A are sliding in A but not in B, and the contacts in 
B are sliding in B but not in A. Let S contain contacts 
that are sliding in both states. We want to show that 

sgn = sgn for at least one a e A n B. (76) 

We will begin by assuming the contrary: 

sgn ^ sgn for all a e A n B, (77) 

and show that this leads to a contradiction. 

The independent load condition Eq. H58|l implies 

Q(v^ - v^, S) = ^ F^^V^ - Y: F^^V^. (78) 

Now let us assume stability for a third state C. The 
following contacts shall be sliding in C: 

• All contacts in §, who are sliding in both states A 
and B. 
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• Some contacts that are sliding in A but not in B. 

Let A' C A denote these contacts. 

• Some contacts that are shding in B but not in A. 
Let B' c B denote these contacts. 

Table Hvl summarizes this information. 

Let C = A U B and C = A' U B'. Note that C can be 
any subset of C. The stability condition for state C is 

Q(v^ - v^,§UC') > 0, 
Q(v-4_v^,§)+^F„^^(V;^-V;^) > 0, (79) 

aec 

Combining this with Eq. H78|l yields 

pAByB _^ pAByA 

E ''^^^ - E P^^'^cf > (80) 

aGB' ael\B' 

Note that A' and B' arc arbitrary, so there is a large 
number of such relations. 

To write these relations in a more compact form, we 
define 



pAByA 
_pAByA 



for a £ 
for a G 



with an analogous definition for 0^. Now the relations 
Eq. 180|l can be written 

E^"+ E ^">0- (82) 

aeC aeC\C' 

And the hypothesis Eq. H77|l becomes 

sgn0;^ ^ sgn0f for aU a e C, (83) 

To show that Eq. (|83|l leads to a contradiction, it suffices 
to show 

E + E 0a > for 1 < * < n. (84) 

where n is the number of elements in C, and Ci is a 
subset of C that contains exactly i elements, and is 
any subset of C^. 

The case i = 1 of Eq. H84() contradicts the hypothesis 
Eq. (113 • To see this, let Ci = {a}. Choosing = in 
Eq. H84|) leads to (f>-^ > 0, and choosing C'l = {a} leads to 
> 0. Thus 0„ and 0f do not have opposite signs, as 
assumed in Eq. (|83|) . This means that Eq. H77() is false, 
and Eq. (|76|l must be true. Eq. (|76|) is equivalent to the 
status change theorem. 

We now show Eq. (|84|l by induction, beginning with 
i = n and proceeding to « = 1. The case i = n is trivial, 
since C„ = C. In this case, Eqs. 182|l and (|84|l identical. 

Now let show that if Eq. (jHH) holds for i + 1, then 
it holds for i also. Choose a contact (3 such that /3 ^ 



Cj, but (3 e C. By the hypothesis, Eq. holds for 
Ci+i = Ci U {/?}. Next, we make two different choices 
for C-_|_;^ and apply Eq. First we choose C'^_^_i — C- 

and obtain 



E'^^+ E ^i>-^l 

and next we choose C'^_^_l = C ■ U {/?}: 



E^"+ E 



(85) 



(86) 



qGC' ael 



Note the parallel between these conditions and Eqs. H67() 
through (O. By Eq. either (p^ or t/)^ is negative. 
Therefore, the only way for both of these inequalities to 
hold is if the sums are positive: 



E'^f+ E '^^>o 



(87) 



This completes the induction step, and thus the proof. 



V. DISCUSSION AND CONCLUSION 
A. How restrictive are the assumptions? 

This work considered only circular particles. The con- 
clusions are probably not modified if other shapes are 
considered. The particle shape most strongly affects the 
geometric stiffness, which is neglected because of the 
quasi-rigid assumption. The moment of inertia plays only 
a small role, because it is eliminated in the quasi-static 
approximation. To accommodate particles of different 
shapes, the particle radius r in Eq. H13|l must depend on 
the contact, and the torque may also depend on the nor- 
mal force. This requires modifying the matrix c. But 
none of these should alter our considerations of mech- 
anism indeterminacy, nor alter the two premises of the 
proof used to resolve contact status indeterminacy. 

Another assumption that seems quite restrictive is the 
use of the linear force law in deriving the stiffness matrix. 
This is no restriction, because this paper revolves around 
the question of what occurs at one point in time, when 
the system must adjust the status of the contacts. There- 
fore, one could always linearize the force law around the 
positions of the particles. 

A second assumption is that all possible states must be 
stable. This assumption is reasonable because if there is 
an unstable state, the packing will probably collapse, ren- 
dering the question of uniqueness irrelevant. The state 
that is most likely to be unstable is the one where all 
Mc critical contacts are sliding. This is so, because the 
contributions of sliding contacts to the quadratic form 
are on the average negative 26]. Furthermore, the only 
way to obtain instability is for negative contributions of 
the sliding contacts to the quadratic form outweigh the 
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positive contributions of the non-sliding contacts [see dis- 
cussion of Eq. (|49|l ] . Now, the packing is always close to 
this state, because when a contact becomes non-sliding, 
it leaves the boundary of the Coulomb cone and becomes 
non-critical. Therefore, when the assumption of stability 
for all possible states is violated, we expect the packing 
to yield. 

Furthermore, if there are vibrations in the system, they 
will be governed by Eq. (|38|l . These vibrations will cause 
the relative motion at each contact to fluctuate. At crit- 
ical contacts, the contact status will therefore switch be- 
tween non-sliding and sliding. Thus the packing will sam- 
ple many different possible states, and if it finds an un- 
stable one, it may collapse. Thus all states must be sta- 
ble in order to guarantee that the vibrations will damp 
out, meaning that this is a necessary condition to obtain 
quasi-static balance. 

Finally, let us remark that theory presented in Sec. If HI 
requires some modifications to deal with opening and 
closing contacts. To account for a contact that opens, 
it is necessary to introduce the status "open" , and allow 
Dn to become positive. Another problem is presented 
by contacts that are initially "open" but may later close. 
In the theory, the particle deformations are assumed to 
be infinitely small, so that no two particles separated 
by a finite distance will ever touch. This may lead to 
the omission of important effects when the particle sep- 
arations are very small, such as in a regular packing of 
almost monodisperse spheres. 

B. Implications of the result 

This work suggests that the stiffness matrix, together 
with the principle of minimum sliding form a complete 
description of quasi-static granular material. Since the 
globally consistent state always exists and is unique, 
there is no need to appeal to other processes that have 
been left out of the model to decide between various pos- 
sible states. Instead of rather brutally setting the particle 
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displacements to zero, as is done in various "stress-only" 
approaches to granular matter, one should consider tak- 
ing the "quasi-rigid" limit, where the stiffness of the par- 
ticles diverges, and the displacements become infinitesi- 
mally small, but are not set to zero. Taking this limit 
leads to the stiffness matrix approach discussed in this 
paper. This work supports the conjecture that the quasi- 
rigid limit preserves all the necessary physics needed to 
describe the quasi-static behavior. 

Of course, there remain some open questions. For ex- 
ample, the question of opening and closing contacts has 
not been dealt with. When a contact reaches the apex 
of the Coulomb cone, it can then go into four different 
states. It can become non-sliding, and move into the in- 
terior of the Coulomb cone. Then, there are two distinct 
sliding states - each one corresponding to a different side 
of the Coulomb cone. Finally, the contact can open. Are 
these four states mutually exclusive, as we have shown to 
be the case for the sliding and non-sliding states in this 
paper? 

This work should also encourage the use of numeri- 
cal methods based on the stiffness matrix. The problem 
faced by these methods is, of course, to find the globally 
consistent contact status. This work shows that such a 
state always exists, and is unique. Thus any way to find 
the state is acceptable. Furthermore, perhaps it is pos- 
sible to use the results of this paper to design intelligent 
strategies for finding the globally consistent state. 
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